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ABSTRACT 


A new algorithm is presented using a logarithmic barrier 
function decomposition for the solution of the large-scale 
multicommodity network flow problem. Placing the com- 
plicating joint capacity constraints of the multicommodity 
network flow problem into a logarithmic barrier term of the 
objective function creates a nonlinear mathematical 
program with linear network flow constraints. Using the 
technique of restricted simplicial decomposition, we 
generate a sequence of extreme points by solving 
independent pure network problems for each commodity in a 
linear subproblem and optimize a nonlinear master problem 
over the convex hull of a fixed number of retained extreme 
points and the previous master problem solution. 
Computational results on a network with 3,300 nodes and 


10,400 arcs are reported for four, ten and 100 commodities. 
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Use EN ERODUCTTON 


Multicommodity network flow problems emerge when several 
distinct commodities flow through a common capacitated 
network and share one or more arcs that are subject to joint 
capacity constraints. The objective is usually to find the 
minimum cost flow given demands and supplies of the com- 
modities. Problems of this type are also referred to as 
"multicommodity capacitated transshipment problems" and 
"multicommodity flow problems" (MCFPs). 

The problem of optimizing MCFPs arises frequently in 
logistic systems. As long as only a single commodity is 
involved, even large-scale problems can now be solved 
routinely by specialized network codes that exploit the pure 
network structure of the problem (e.g., GNET by Bradley, 
Brown, and Graves [Ref. 1]). Such solvers are not directly 
applicable to the general MCFP, and general linear program 
solvers are usually inapplicable as a result of the large 
constraint matrix encountered with typical MCFPs. 

Because of the importance of the MCFP, much effort has 


been devoted to finding efficient, specialized solution 


techniques. Earlier surveys are given by Kennington [Ref. 2] 
and Assad [Ref. 3]. New and effective decomposition methods 
were recently developed by Staniec [ Ref.4] and show 


encouraging results. This paper extends that research by 


deriving and implementing a decomposition for the MCFP based 


on a logarithmic barrier function. 


A. 


the 


STATEMENT OF THE PROBLEM 

In order to formulate the MCFP as a mathematical program 
following notation is used 

G = (I,J) is a network with set of nodes I and set 


of arcs J. 


P is the set of commodities (products) flowing on G. 
= e I is a node of G. 
j e J is an arc of G. 


p € P is a commodity flowing through G. 


No is an |I| x |J| node-arc incidence matrix for each 
product (Ny =N>z= ... = Nip) - 
N is an |I|°|P| x [J]*|P] matrix with matrices No 


along the diagonal, Os elsewhere. 


A is a (Jie a matrix. (lo ieee. 


C = (Cq1---+C¡py) is a vector of arc costs, length 
PA 

X = (Xqr or X py) is a vector of arc flows, length 
e 


by is a vector of joint capacities with length |J]|. 
Dx is the vector (Dir - + -7D3) with length |J|°|P]. 
by 1s the vector of Supplies and demands for each 


commodity with length |J|°|P|. 


Then the MCFP may be stated as follows 


(P) min CX (duals) (ae) 
S.l Ax <S Dy (uy) (2) 

Nx = D3 (us) (3) 

Osx Say (u3) (4) 


For the constraint sets (3) and (4) the abbreviated notation 
x €e F will be used. 

The stipulation that binds the flow of several com- 
modities to joint capacities is given by (2) and constitutes 
the set of "complicating constraints". Without those con- 
straints, the problem would reduce to independent, bounded, 
single-commodity flow problems. The duals u,; corresponding 
to (2) are nonpositive. The constraints (4) are redundant, 
but they ensure that x is bounded when the constraints in (2) 


are relaxed. 


B. SOLUTION METHODOLOGY 

The two basic approaches to solving the MCFP (other than 
a Standard primal linear programming method) can be charac- 
terized as either decomposition or partitioning techniques. 
The latter employ a special basis factorization within a 
Simplex algorithm in such a way that portions of each 
generated basis maintain characteristics of the pure network 
flow problem. Those method are not investigated here. For 


further detail see Kennington and Helgason [Ref. 5]. 


Decomposition methods solve MCFP by using a master 
problem that coordinates the solution of single network flow 
subproblems. These methods are attractive, since they may 
require the internal storage of only one commodity at a time. 
This approach can further be divided into resource-directive 
and price-directive algorithms. Both will be stated here for 
later reference. 

The resource-directive method uses a master problem that 


distributes improving capacity allocations to the individual 


commodity subproblems. For this purpose MCFP may be written 
as 
(P’) ee min Cx (5) 
Se ay b; (6) 
x-y<so0 (7) 
Nx = b, (8) 
DOSIS ba (9) 
SES Dig (10) 


Any vector y satisfying constraints (6) and (9) may be 
interpreted as a "capacity allocation" which apportions the 
capacity of an arc across the individual commodities. The 
inner minimization for a fixed y amounts to a restriction of 
(P). If its solution is feasible in (P), it yields an upper 


bound on the optimal solution. 


thes sibproblemss forsas particular. allocationmvector y are 


of the form 


(SP1 (Y)) V(y) = min cx (11) 
Sete Nx = Do (12) 
0OSxS<S y (E) 


which is a set of  single-commodity minimum cost flow 
problems, solved independently for each commodity. 


The master problem than becomes 


(MP1) min V(y) (14) 
sker Ay < b} GES) 
OS Sere (16) 


The objective function in (14) is piecewise linear and 
convex. In theory, the master problem can be solved by 
subgradient optimization or by a cutting plane algorithm. 
Penalty and barrier function decomposition are examples 
of price-directive decomposition. Before describing them, it 
is worthwhile looking at the Lagrangian relaxation of (P). 
If the joint capacity constraints are placed into the 
objective function with multipliers uj S 0, the Lagrangian 


dual of (P) is 


(LR) max min EU a uy (by - Ax) (17) 
u, 380 x20 


St. x € F. 


This problem corresponds to solving max LR(u,) where 


uy $0 
LR(u,) 1s defined by 
(LR (u, ) ) min cx - uy (Ax - b) 
xEF 
= min (c - uy A) x + u,b,- (18) 
xeF 


Note that the evaluation of LR (u4) requires the solution of 
|P| independent single commodity problems. Furthermore, for 
any fixed uz = 0; LR (u, ) ylelds a lower bound on (P) and this 
bound will Þe tight If u is optimal in (P): 

LR (u; ) = a 

LR (uy) is a piecewise linear and concave function that 
can be optimized, in theory, by subgradient optimization (for 
example, see Fisher [Ref. 6], Goffin [Ref. 7] or Sandi [Ref. 
8]), or by a cutting plane method ( see Kelly [Ref. 9]). 
Optimization of LR(u,) by a cutting plane algorithm is 
essentially equivalent to solving the MCFP by Dantzig-Wolfe 
decomposition (see Staniec, [Ref. 4]). 

However, even if (18) is solved optimally, it is possible 
that LR (uy) = rl for Uy £ oe Furthermore, we may not be 
able to obtain a corresponding primal solution that is 
feasible to problem (P). 

Solving LR(u,) for any u, generally allows some 
constraints to be violated with penalty -uy (Ax - by). Other 
penalty functions are possible which will, at least in a 


limiting sense, yield optimal primal solutions. The 


objective function in (P) can be transformed into a nonlinear 
auxiliary function Q(h,x) that typically includes a polynom- 


lal penalty term of the form 


h 
cx + — |] (Ax - 2 E 


EEN "On X) E z 


xEF 


5) 


ex + q(h,x) (20) 


where || ° lilies indicates the ¿th 


vector max(0,Ax - b), i.e. the vector of violations of the 


norm, and (Ax - bi), is the 


constraint set (2). The value of h constitutes a positive, 
increasing sequence of penalty parameters. Furthermore, 
t > 1 is required for this method to ensure convergence. 

The penalty function has some attractive properties, as 
given in Luenberger [ Ref. 10] : Let { nk } and { xk } be 


sequences of penalty parameters and optimal solutions to 


(PP (hX)), respectively, with nxt 2 nk and n? > 0. Then 
Q (hË, x5) < o(nX+1 ¿X+1) (21) 

q (hf, x£) 2 an? ti) and (22) 

cxk < cxktl | (23) 


The algorithm converges to the optimal solution. Thus, 


* * 
lim Q (hË, x5) = cx , and xk > X 
h*—=0 


Since a feasible solution is usually not obtained until final 
convergence occurs, the penalty approach is classified as an 


exterior point algorithm. 


In order to obtain intermediate feasible solutions a 
suitable scaling or projection procedure may be used. If we 
model the MCFP to include additional bypass arcs which 
satisfy any undeliverable demand at a high cost from a 
supersource, we can assume that any allocation of capacity 
satisfying (6) and (9) is also feasible in F. One pos- 
sibility to obtain such capacity allocations y from a given 


vector x e F is given by defining y = CA (x) as 


A A å A ine > 

3b, 4 / (AX) if (Ax-b,) 5 2 0 
Ypj e 
A E E EA DE 0. 


Then, a solution of (SP1 (y)) as defined in (24) yields a 
valid upper bound at the current iterate x. This allocation 
procedure was successfully applied by Staniec [Ref. 4]. 

The idea behind the general barrier function approach is 
just the opposite of the exterior point penalty methods. 
Starting with a feasible solution which lies within the 
interior of the constraint region, a modified objective 
function establishes a barrier against leaving the feasible 
region. For the MCFP, the auxiliary function with respect to 
the joint capacity constraints then becomes 


(BP (h) ) min  B(hN,x) = cx f BED X). (25) 
xeF 


An ideal barrier term b(x) would take the value zero for 
all interior points and infinity at the boundary. A sequen- 
tial barrier function b(h,x) omits this discontinuity and may 
take the following forms 


ae (26) 


Jo Oar) 
or 


b(h, x) ZO 2. ln (b4 = AX) jv CLT 


h component of the vector 


where (by = Ax) y denotes the qt 
(by SSB 
The first expression (26) is called the "inverse barrier 


function" and the second term (27) the "logarithmic barrier 


function" (see, for example, Fiacco and McCormick [Ref. 11] 
or Bazaraa and Shetty [Ref. 12)). Both functions approach 
the ideal barrier function b(x) as h > 0. Any barrier 


algorithm requires the existence of an interior region, i.e. 
it does not work for equality constraints. 

The logarithmic barrier function was first proposed by 
Frisch 1955 [Ref. 13]. It was then derived together with the 
inverse barrier function from the Kuhn-Tucker conditions for 
optimality by Fiacco and McCormick (Ref. 11]. The logarith- 
mic barrier function has obtained recent attention in linear 
programming due to its fundamental properties. Megiddo [Ref. 
14] investigates the properties of a weighted logarithmic 
barrier function for general linear programs that places the 
nonnegativity constraints on x into the barrier term while 


requiring strictly positive values for x. This approach 


leads to a smooth path through the interior of the constraint 
region towards the optimal solution and theoretically 
validates its use in linear programming. Gill, Murray, 
Saunders, Tomlin and Wright [Ref. 15] established a close 
relationship to the projective linear programming algorithms 
initiated by Karmarkar [Ref. 16]. 

The transformation of the objective function into a 
penalty or barrier function creates a nonlinear programming 
problem which requires an efficient solution method to make 
such transformation profitable. Staniec [Ref.4] successful- 
ly applied the method of restricted simplicial decomposition 
(RSD) in order to solve the penalty function decomposition. 
The idea was developed by Hearn, Lawphongpanich and Ventura 
(Ref. 17]. RSD solves any pseudoconvex optimization problem 
with linear constraints by generating extreme points in 
linear subproblems while a master problem optimizes the 
Original objective function over a simplex derived from a 
fixed number of retained extreme points plus the last 
iterate, i.e. the last solution to the master problem. The 
implementation of this method for the MCFP will be described 
later in more detail since it proves useful for the barrier 


decomposition as well. 


C. TEST PROBLEMS 
A real-world large-scale MCFP should be used in order to 


examine the efficiency of the proposed algorithms. Staniec 


10 


developed an appropriate problem that describes the 
transshipment of conventional ammunition from production and 
storage locations to overseas debarkation points and theatre- 
Of-war locations via capacitated road, rail, sea and air 
eransportation links; The product demands are time-phased 
and the objective is to minimize the weighted deviation from 
on-time deliveries. The network contains backlogging arcs 
with graduated penalties and bypass arcs that satisfy 
undeliverable demand at high cost. For more details see 
Staniec [Ref. 4 and 18]. The same network is used to test 
and compare the algorithms for four, ten, and 100 commodity 
problems. The underlying network contains 3,300 nodes and 
10,400 arcs of which about 10% are subject to non-redundant 


capacity constraints. 
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II. THE CONCEPT OR T BARRIER FUNCTION TDPDECOMPOSTILION 


The original MCFP (P) constitutes a linear program. It 
is only the size of the problem that makes a primal simplex 
algorithm unattractive and computationally expensive. 
Penalty and barrier functions were initially designed for 
nonlinear programming problems where they proved useful in 
converting a constrained nonlinear optimization problem into 
an unconstrained nonlinear problem. For the special case of 
the MCFP, penalty and barrier functions provide good 
decomposition tools. 

The barrier function decomposition for the MCFP retains 
the basic decomposition idea of the penalty function 
approach. The overlapping joint capacity constraints in (2) 
are placed into the barrier term of the objective function 
in order to enable the successive solutions of independent 
Single commodity problems in a sequence of subproblems. The 
selection of the logarithmic or inverse barrier function is 


not arbitrary but has an interesting theoretical derivation. 


A. THE DERIVATION OF THE BARRIER FUNCTION 

Fiacco and McCormick derive the barrier function from 
the Kuhn-Tucker sufficiency conditions for constrained 
minima [{Ref. 11]. A good and detailed analysis, including 


implementations and numerical results, is further given by 


12 


Wright [Ref. 19). The derivation shows that the use of a 
logarithmic barrier function is not arbitrary since it has a 
very natural origin. 

For the purpose of this analysis it is convenient to 


restate the problem (P) in the following form 


(pes) Min f(x) = cx (28) 
XEF 
Slate digg = by > Ax = 0 (29) 
where the objective function (28) has gradient Yf(x) = c and 
each constraint in (29) has gradient Vg 4 (x) = - ad, the 


negative of the ¿th row in A. We associate the dual vari- 
ables M4 with the constraints in (29) and note that in 
comparison with the duals u} of problem (P), M} now takes 
the opposite sign : Uy = - Hy- 

Following the derivation of Fiacco and McCormick, we 
assume that there exist points in the neighborhood of the 
optimal solution to (P’’) such that strict inequality holds 
for the constraints in (29) i.e., g(x) > 0. Furthermore, we 
allow a perturbation of magnitude h in the Kuhn-Tucker 
Suc ieney Conditions fOr Optimality. At some point 
[x (nh), uy (h) ] near the optimum (x",u, ) the following 
conditions have to hold for h small and for all j € J: 

(by - Ax) y > 0 (primal feasibility) (30) 
Hij (by - Ax). =h > 0 (31) 


J 
(perturbed complementary slackness) 
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IV 


0 (32) 
(dual feasiblity) 
Ox: (33) 


Pas 


= , 0 
C 25 His (-a-) 


Rewriting (31) as His =h / br Ax) y and substituting in 


(33) yields 
= 1 - pe 
c - h Xs [(by - Ax) 5] (-aJ) = 0. (34) 
Using the notation gj = (by - Ax) 5 20- forall je 17 


equation (34) is of the general form 


V9 5 [x (h) ] 


il 
© 


YB(h,x) = Vf(x(h)) - h 2, 


j (35) 


gj[x(h)] 
and simply means that the gradient of the objective function 


B(h,x) = f(x) - h 2. In (be) (36) 


J 
vanishes at x(h). This is the logarithmic barrier function! 
Note that no constraint qualifications are necessary since 
all constraints in the MCFP are linear. 

Fiacco and McCormick show that the second order 
sufficiency conditions are also satisfied for B(h,x) at 
(x (h) , uy (h) ] near the optimum xo and prove the existence of 
x(h) satisfying these conditions. It is also shown that the 
Hessian matrix is positive definite for small h. 

The same authors obtain the inverse barrier function 
from a modification in the derivation above that enforces 
nonnegativity in (32) by introducing a variable T such that 
2 


tT = Hij: The logarithmic barrier function seems to be the 


more fundamental approach. 
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BBS PROPERTIES OF THE LOGARITHMIC BARRIER FUNCTION 

The basic properties of the barrier function are again 
given by Fiacco and McCormick [Ref. 11] as well as by Wright 
(Ref. 19] and can be stated as follows for the MCFP 
For a decreasing sequence (nå) and associated minima {xk} 


the following conditions hold 


B(hX,xX) <s pB(nT1 xxl, en) 


for sufficientiø small nk and bounded (by - AX) 4, and 


cx < cxkrl (38) 

; - Eer. = Pel 

2, 1n (by Ax) y < 2. 1n (by Ax) y p (39) 

. k k 

lim h* 2, In(b, - Ax).* => 0, and (40) 
ca, J : J 
Lim da e (41) 
h*=>0 


The following small example will illustrate these 
properties. The problem consists of a nonlinear objective 
function and a single constraint : Min 2x? s.t. x 2 1 and 
has the optimal solution x = 1. The barrier function 


B(h,x) = 2 x? - h ln(x-1) has a closed form solution 


Os + 40.25 4M0N25 hy) °°? 
The solutions are shown in Figure 1 for linearly decreasing 
values of h from 20 to 0.1. The change in x as a function 


of h and the corresponding logarithmic term are depicted in 


Figurer? 


L5 


FCXCH)) 
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3(X,H) 


X(H) 
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16 


Eewanccm tts boundary. => Mé is Smooth and can be well 


approximated by a straight line until B(h,x) reaches its 


maximum. From that point on the rate of change in x(h) 
increases slightly. Similar results are obtained for 
different objective functions. Returning to the multiple 


constraint barrier function for the MCFP, the maximum of 
B(h,x) is obtained for a value h’ which satisfies the sta- 
wonarity condition 2¿1n (by - Ax (h'’)) j= Os This scan only 
occur if some components of (oie en) are less than or 
equal to one, i.e., some constraints are almost tight. 
Thus, we will not be able to observe the property in (37) 
until the finai stage of the algorithm when relatively small 
values of h are attained. 

The properties (38) and (41) are the most important ones 
for practical purposes. Starting with a feasible (interior) 
point, the algorithm produces a nonincreasing sequence of 
objective values which converges to the optimal value of (P) 
and provides intermediate feasible solutions along its path. 
This path of x as a function of h describes a trajectory 
that has already been studied by Fiacco and McCormick [Ref. 
11] as well as by Wright [Ref. 19]. The existence of this 
trajectory could be utilized in an extrapolation technique 
predicting the next iterate or even x, Its implementation 


for the large-scale MCFP is not investigated here. However, 
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we will be able to get a good feasible solution after only a 
few iterations without using an extrapolation technique. 

One final property is the robustness of the barrier 
function method: Mifflin (Ref. 20] showed that it is 
sufficient to solve B(h,x) at each iteration only 
approximately, within a predetermined tolerance, while still 


achieving convergence (at a lower rate). 


C. COMPUTATIONAL CONSIDERATIONS 

Any barrier function technique requires an interior 
Starting point. It may not be easy to find a starting point 
and the performance of an algorithm is influenced by the 
quality of this initial solution. We utilize the capacity 
allocation mechanism y = CA (x) for the generation of an 
initial starting point. 

The transformation of (P) into a nonlinear programming 
problem requires an effective NLP solution methodology. If 
a line search is part of this method, any discrete-step line 
search procedure along an improving direction may cause the 
evaluation of B(h,x) at one or more infeasible points. This 
requires extra precautions during the implementation, 
resulting in additional computation time. 

Ryan [Ref. 21] points out that for small values of h the 
auxiliary function B(h,x) becomes very "steep valleyed" and 
the gradient VB(h,x) can take large values in a small 


neighborhood of x(h). Therefore, a termination criterion 
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based on the magnitude of the gradient alone may be critical 
as h becomes small; termination criteria based on the 
difference in successive objective values may lead to 
premature termination. We will consider both criteria and 
take the risk of not solving the problem optimally at each 
step. 

Another well-known difficulty oof barrier function 
approaches is the ill-conditioned Hessian matrix of B(h,x) 
for small h (see,for example Bazaraa and Shetty [Ref. 12], 
Wiighte~ Ref. 19] or Ryan (Ref. 21)). This problem emerges 
only in the final stage as h > 0. For large-scale program- 
ming purposes the achievement of a solutiom which defers 
from the optimum only by some small value € ( E€-optimality) 
is normally sufficient. Such a solution may be obtained 
before the ill-conditioning becomes bothersome. 

Another issue that significantly influences the perfor- 
mance of the algorithm is the choice of the initial barrier 
parameter h and its rate of decrease. Lootsma [Ref. 22] 
showed that the absolute difference || x(h) - x Mb Ls TOR 
the order O(h) for the logarithmic barrier function. Ryan 
(Ref. 21] uses this linear relationship to propose a 
generating relation of the form hk = 1017K h? , where 
k = 1,2,3,..., and h® is positive. 

The suitable choice of the initial value nh? is a more 


critical issue, since theoretically h can take any value 


ro 


greater than zero and up to infinity. Intuitively however, 
h% should depend on the cost and constraint structure of the 
problem. One possible choice would be to interpret the 
parameter h as a scaling factor between the cost vector and 
the set of constraints placed into the objective function. 
This leads to the suggestion that we achieve initial balance 


O such that cx? = h? 2; Ln (by ~Ax”) 5. 


at the starting point x 

Instead of analytically deriving hÔ, simply taking an 
multiple 6 (6 greater than one) of the maximum cost value in 
connection with a constant rate of decrease has performed 
well in test problems 


k+1 


h =o0o*c >, h = 2 E foramen 0 e a 


max 
In general, we recommend choosing h? too high rather than 
too low as the algorithm will more quickly adjust to high 


values rather than low values. 
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III. THE LOGARITHMIC BARRIER FUNCTION DECOMPOSITION 


USING RESTRICTED SIMPLICIAL DECOMPOSITION (RSD) 


The basic idea of the barrier function decomposition for 
the MCFEP is to place the coupling constraints (2) into the 
logarithmic term of the barrier function as derived in 
Chapter 2. The resulting formulation (BP (h)) constitutes a 
nonlinear programming problem with linear network flow 
constraints (3) and (4). Using the restricted simplicial 
decomposition technique (RSD), we will decompose the problem 
into a nonlinear master problem and a set of subproblems, 
which require only the solution of |P| independent pure 
network flow problems. The master problem has a reduced 
search space described by a fixed number of retained extreme 
points, which are generated in the subproblems. 

Lower and upper bounds on the optimal solution to (P) 
can be easily established. The analogy with the penalty 
decomposition is interesting and worth pursuing. The 
solution method RSD used in the penalty function decomposi- 
tion also proves to be effective for (BP(h)) and is descri- 


beds first 


A. RESTRICTED SIMPLICIAL DECOMPOSITION (RSD) 
The basic difference between a linear programming 


problem and a nonlinear programming problem with linear 


2l 


constraints is the fact that the optimal solucrones,... 
normally not be an extreme point of the constraint region. 

The familiar Frank-Wolfe algorithm [Ref. 23] takes 
advantage of the specialized constraints by generating an 
extreme point solution of the original problem in a linear 
subproblem whose objective function is the linearization of 
the original objective function at the current iterate 
(given an initial solution). A master problem provides a 
new iterate via a simple line search between the previous 
iterate and the new extreme point. The main disadvantage of 
this decomposition algorithm is its susceptibility to slow 
convergence, especially when the line search direction 
becomes orthogonal to the gradient of the objective function 
(e.g., see Wolfe [Ref. 24]). 

The method of simplicial decomposition is due to 
Geoffrion [Ref. 25], von Hohenbalken [Ref. 26] and Holloway 
Ref. 27]: A nonlinear master problem replaces the line 
search of the Frank-Wolfe method by extending the optimiza- 
tion to the convex hull of all extreme points generated in 
the linear subproblems. This solution method relies on an 
effective solution of the nonlinear master problem. 
Although the master problems have only simple convexity 
constraints, the increasing number of extreme points makes 
the implementation of this technique unattractive for the 


solution of large-scale programming problems. 
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The disadvantage of simplicial decomposition led to the 
idea of restricted simplicial decomposition (RSD) as develo- 
ped by Hearn, Lawphongpanich and Ventura [Ref. 17]. RSD 
limits the size of the master problem by fixing the maximum 
number r of retained extreme points. The master problem 
optimizes over the simplex of these extreme points and the 
iterate obtained from the previous master solution. Any 
newly generated extreme point replaces the old extreme point 
with minimum weight in the expression of the current iterate 
as a convex combination of the retained extreme points and 
the prior iterate. After solving the new master problem, 
all extreme points with zero weight can be discarded. 

If r is set to its minimum value r = 1, RSD specializes 
to the algorithm of Frank-Wolfe. For the maximum value of r 
(the finite number of extreme points) the method represents 
Simplicial decomposition. The solution to the master 
problem becomes harder as r is increased, but is rewarded by 
Significant improvements in the convergence rate to the 
optimal solution (see Hearn, et al., [Ref. 28]). 

The decomposition of (P’’) is formed in the following 
manner: Let XË denote a matrix in the kth iteration of the 
master problem whose columns are a set of r extreme points 
from F, let xk-l be the solution of the previous master 


k=1, i 


problem, and let xX = xk y {x The master problem in 
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terms of the weigħts w at iteration k been 


for a fixed h, 


(Mp2K) min B(h,X*w) = cxkw - h E, ln(by - AXKw) j (42) 
s.t. DL. (43) 
w>0, (44) 


which has solution wk in terms of w and solution xk = xkyk 


in terms of x. 
The subsequent subproblem optimizes the linear ap- 


proximation of B(h, z) E xX over F, which is equivalent to 


(sp2K) min 9B(h,x*)x . (45) 
XEF 
It is convenient to introduce the notation (b, - Ax)t to 


represent the vector (b, - Ax) whose components are taken to 


ER 


the t power. Then, (sp2X) can be written as 


(sp2K) min (c - h{(b, - AR) A)x . (46) 
XEF 


Using the relationship in (31), we observe that we obtain an 
estimate Ha of the optimal dual variables as = h (by -Ax*) 7? 
at each iteration. Substituting into (46) yields a sub- 


problem as 


(SP2*) min (c + [l,A)x. (47) 
XEF 


This subproblem resembles a standard Dantzig-Wolfe decom- 
position subproblem as it would be used in a dynamic column 
generation approach to problem (P) with dual estimates 


A 
u = - H3- This issue is not discussed here; see Staniec 


(Ref. 4]. 
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DUER EOMENSRSTOCERSErFUS ture Fof FEhe constraint matrix N, 
(one) permits independent solutions for each commodity. 
However, each solution to (sp2X) is not necessarily a 
feasible flow in (P) and of course, is not an interior 
pornt: AS a starting point, the barriér function approach 
requires an initial interior point, which may be obtained as 
follows. First we solve a subproblem with h = 0, and 
yielding solution x0. Then, using the capacity allocation 
mechanism Y = CA (x0), we further reduce the capacity by 
setting y" > bær y fon O b <el: Finally, we solve 
subproblem (SP1(Y")) to get an initial interior point 
saint zon to (F). 

Starting with this solution in the master problem we 
simply limit our search to that part of the convex hull of 
extreme points plus the current iterate which provides an 
interior point solution as next iterate. This will be 
discussed in connection with the solution of the nonlinear 


master problem. 


B. SOLUTION OF THE MASTER PROBLEM 

The reduced gradient method is a well-known approach to 
solving a nonlinear program with linear constraints and is 
used here for solving the RSD master problem. The method 
was first proposed by Wolfe [Ref. 28] and modified by 
McCormick [Ref. 29]. The basic idea is the partition of the 


variables x into m basic variables Xp and n nonbasic 
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variables x, as done in the simplex algorithm of linear 
programming. This induces a partition of the constraint 
matrix A into parts B and C, where B is assumed nonsingular. 


The NLP then takes the form 


min f(x) = £ (Xp, Xy) (48) 
s.t. Ax = BXp + CXy = b (49) 
XgrXy 2 0. (50) 
The variables Xy are regarded as independent variables 


whereas x, are dependent variables completely determined by 
Xy and equations (49). Consequently, the objective function 
can be considered to be a function of x, only and the 
constraints reduce to the nonnegativity constraints on the 
independent variables and the limitation in their change 
that provides nonnegative basic variables. This fact allows 
the application of a modified steepest descent method 
accounting only for the nonnegativity constraints. The 


k 


reduced gradient r^ at iteration k is of dimension n-m and 


computed as 


k 


D k 
r` = (rB 


E Pac £ (Br XN) z Ix f (Br XN) Bic. (51) 


To find an improving direction d such that v£(x%)d < 0, we 


select at each iteration k 


dyj — for Xy nonbasic (52) 
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and 
dai = - (B` Cdy? i for. basic (53) 

A new direction needs to be computed as soon as a nonbasic 
variable attains its zero level. If a basic variable 
becomes zero, the partition must be modified. Also, this 
method requires a nondegenerate solution at each iteration. 
A more detailed description can be found in Luenberger 
[Ref. 10] or Bazaraa and Shetty [Ref. 12]. 

The use of the reduced gradient method for the solution 
of the master problem (mp2X) results in some nice simplif- 
ications due to the presence of only a single convexity con- 
straint. There is only one basic variable w, to be selected 


and the reduced gradient becomes 


OS acta (bh ES (Ax: - Axį) (54) 
for the nonbasic variables Wj: Computing the direction 
component for the basic variable w; reduces to 

A A IS A 
d; B Cd, Xd, (55) 
and for the nonbasic variables W 
= 8 SL ie D Ss 0 
= 56 
dj (56) 


Zwi ; o> 4 
War, if rj 0 
The line search in the direction d is limited by the non- 


negativity constraints. Thus a maximum steplength QU” is 
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computed as 
ao’ = min {-w,/d, | dy S10; tor all varraciles Wg}. (57) 
Any move in the direction d of size @ must also satisfy 
by - [AX(w + ad)] > 0 for allm Ums gS OTI 
or equivalently 
Dy - AXw > QAXd for all 0S A SA”. (99) 
Since by - AXw > 0, this holds for all arcs 1 with 
(AXd) 1 < 0. EE (AXd) y > 0 for some arcs l, we perform a 
ratio test to select 
ao’ < min ( (by - AXw)y (AXd), | (AXd), > 0) (60) 
and set 
Aan min Ol’ , OL). (61) 
The master problem solution is summarized as follows: 


Step O : Set w such that lw =1, Wh 2 0. 


Select a basic variable: use the largest w,. 


Step 1 : Compute the direction d as determined in equations 
(54) through (536). tit d= Uy seop. 
Step 2 : Solve min B[h,X(wtad) ] SE. TONES Q S Cae ina 


line search, yielding @,. Let w <— w + and and 
go to Step Ik 
The convergence of this method may not be satisfactory 
since it represents a greedy descent method in the reduced 
space of the nonbasic variables. Reklaitis et al. [Ref. 30] 
suggest a convergence acceleration technique by using a more 


effective unconstrained search method as long as the basic 
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variables do not change. This is especially relevant in our 
case where the prior iterate frequently takes the largest 
weight and remains the basic variable. Therefore we apply a 
conjugate gradient method, that has proven its effectiveness 
in unconstrained optimization and is easily implemented. 
The modified direction for the first iteration and any 
iteration following a basis change becomes 

= rå me Wj > 0 or ab < 00 

yz = (62) 

0 otherwise, 

and for all other iterations 
ió 


k k k-1 
== — e AAA AAA A A 
d N E z a I dy (63) 


The direction for the basis variable is the same as in the 
standard reduced gradient method. The use of the conjugate 
gradient provides significant improvements. In a typical 
master problem with only four extreme points, the reduced 
gradient method used 500 iterations to obtain a solution 
that was still 13% worse than a solution obtained with the 
conjugate gradient method in 28 iterations. 

Each line search requires frequent evaluations of the 
objective function and consumes a fair amount of computation 
time. We can improve this by modifying the objective 
function to include only those arcs in the barrier term that 


are potentially able to violate the joint capacity con- 
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straints at the current iteration. These arcs are easily 
identified and recorded in a set Jy containing all arcs that 
ever had a violation in any of the extreme points generated 
in the subproblems. The arcs j é Jy cannot violate the 
joint capacity constraints in a solution to the master 
problem that is a convex combination of the extreme points 
and the feasible prior iterate. Thus we establish the 
barrier only on a reduced subset of the joint capacity 
constraints. This set is updated at each solution of a 
subproblem in case a new arc experiences a capacity viola- 
tiong This procedure amounts to a modified barrier func- 
tion. However, in a finite number of iterations, the number 
of arcs in Tyr will take a fixed value | Jy | < |J| and no more 
arcs will be added to Ją. From this point on we are back to 
a pure barrier function as derived in Chapter 2 and conver- 


gence of the algorithm is preserved. 


C. LOWER AND UPPER BOUNDS 

Since we will rarely be able to find the optimal 
solution to (P) within a reasonable number of iterations, we 
need to establish bounds on the optimal solution. 

Lower bounds can be derived from the Lagrangian dual 
problem 


(LR (uy) ) min (c - u¡ A) x + u,b, , which provides a lower 
XEF 


bound on (P) for any fixed u] S 0. 
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A 


Using the dual estimates uy = - Hy, E Ze 


u =~ h (By 5 Ax) 1 for each solution x provided by the 


master problem, we obtain 


(LR (44)) min a 4, A) x - br (64) 
XEF 
Recalling the subproblem (SP25): min(e r y A) x , we find 
xeF 


by comparison that both objective functions differ only by 
the constant term bb, and yield the same optimal solution 
xÝ. Thus, we obtain a valid lower bound ben by subtract- 
ing the constant term Ub, 

v(u,*) = (¢ + u,*ayx® - tb... (65) 
Furthermore, since in the limit xk > xo and es => gg IE 
must be that U D (Ci ly A)x E Hy by = cx. 

Upper bounds on the optimal solution are generated in 
each master problem, since we restrict the solution to an 
interior, feasible point. Thus, if xk = xkwk solves (MP 2*) , 
the upper bound 

V(x*) = cx* (66) 
is readily available at each master problem solution. Due 
to the convergence property (41) of the barrier function, 
it must be that V (xk) = ose? as h 0. 

However, we are usually able to obtain a better upper 
bound by utilizing the capacity allocation mechanism 


k 


again. Let y = CA (xk) at some iteration k and let x” be an 


Sel 


interior point. Then, forall ees J€ Tyr (Ax - bi); < 0 


A A 


Yo; 2 i for all p,j from equations (24). Then 


the following relationship must hold: 


and 


== NR A AS 
V(x) = cxk 2 mMin Cx 2 Mae mee) (ya 
xEF XEF 
Ke xSy 
Since y 2 xk, 


Thus, solvamnig SP1 (y) yields a feasible solution with value 


Viy) <S xk). 


D. THE ALGORITHM RSDPxXE> 

The algorithm RSD(B) using a barrier function decomposi- 
tion can now be presented. An initial lower bound is 
obtained by solving the problem without the joint capacity 
constraints (2). If this solution is feasible in (P), it is 
optimal. Otherwise we obtain a feasible solution via the 
capacity allocation mechanism, which gives an initial upper 
bound. The algorithm generates a sequence of extreme points 
in the subproblem and an interior point solution in the 
master problem until e-optimality is achieved. As a 
heuristic, we invoke the capacity allocation mechanism at 
every rth iteration of the master problem to improve the 
current upper bound and decrement h at this time, even if 


the master problem is not solved optimally for the cur- 


rent h. 


82 


Notati 


on 
matrix of retained extreme points at iteration k 
current master problem solution 
previous master problem solution 


k 


matrix X” augmented with XM 


optimal solution to subproblem (SP25) 


H (X) convex hull of X 


(x) a capacity allocation based on x 


set of joint capacity arcs which are violated in 


at least one subproblem solution x, 


barrier parameter used in pth 


parameter update 
maximal number of retained extreme points 


stopping criteria for near-optimality. 


AlgorrtnhmeRSsD (B) 


Input 


Output 


Step O 


The network G = {I,J}, Joint capacity vector Dy, 


cost vector c and supply/demand vector b5. 


K and final bounds FV, Vv. 


Best obtained solution x 
(ina tral Zae1on) 
Select hn17 > 0, € 280, r > 1, 0 8b < 1. 


Set k = 0, 1 = 0, 


x? =ø, Jy =O, AS = 0 for all j e J. 
Solve (LR(Q)) : x9 = argmin {cx | x € F }. 
K 
S AQ _ i E A 
SEE NNE KOR = Set Jy ~ (J T (b] AX ) j < 0}. 
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0 


If J,, = 0, stop with x optimal. 


Else, set y = ca(x0), y" = + y and 


solve SP1(y") yielding o 
Set V = exe If (V - v)/V < €, exit with sk ME 
Else, set xl = Ø , Xm = xt = a , k= 1. 
Step 1 : (Solve sp2*) 
Solve xk = argmin {(c + ly FA) x I xe F}, 
x 
k 


z 
Fy 
0 
ES 
© 
po 
Re 
LJ 
Il 


Bior Axo e ae 


HS = 0 for all j € Jy. 
Set Gp Jy U a E 
So OA 


E ER 
TE (PE) e ica << mace 
If (c + pt, “ay (xX - x*) 2 0, x® solves B(ht,x) . 
set htt} =a * nt (0 < asl). 
Set 1 = 1+ 1 and go to Step 2. 
Else 
G) af Ixky < r , x5tl = xE U {x£}. 
(1i) alse ¡xx = r , drop the column of xk which had 


the smallest weight w£ in the convex com- 


i 
k k 


bination forming x^“ and replace it with x^. 


Go to Step 2. 
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Step 2 : (Solve the master problem MP2X) 


set RKHL = xk+l Y {xy}. 
Find xktl = argmin { B(ht, x) | x € H(X5t1) } 


x 


=> Ma where 1 <i < jx*t1) ana 


de 
SOM is the ith column in yktl 
_ ?£k+1 
Set Xm = X 
tf cx®tl 6 UV, set V = cxktl, 


If (V - V) / V< €, exit with x, V, V. 
If k+1 is an integer multiple of r, do a capacity 


KFI 


allocation : Set y = CA(x ) and solve SP1 (Y), 


yielding V(y) = cx 


ie (SV eset x = x,y 


exit with x“, V, V. 
ana aaa a 1. 


Set k = k + 1 and go to Step 1. 


C3 IMPLEMENTATION OF RSD (B) 

The algorithm RSD(B) has been coded in FORTRAN which is 
still an extremely efficient language for mathematical 
programming purposes (see MacLennan [Ref. 31]). A sophis- 
ticated data structure is used for the storage of sparse 
matrices and vectors. Allowing direct communication with 


the X-system solver of Brown and Graves [Ref. 32]. Integer 
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arithmetic is performed to the extent possible, especially 
for the subproblems which use a GNET Solver [Ref. 1] to 
produce rapid solutions. 

A design goal was set to provide a decomposition 
procedure which solves only a single commodity problem at a 
time, and thus operates easily within a modest memory 
region (say, two megabyte virtual storage capacity). Other 
information such as current incumbent, previous iterate, 
prior extreme points, etc., have to be kept on external 
storage devices. This approach leads to considerable 
input/output operations at the expense of computation time, 
but the maximum problem size in terms of number of com- 
modities remains independent of the available virtual 
storage. Also, the resulting algorithm is highly parallel 
by commodities. 

The implementation of the master problem contains 
several parameters such as the number of retained extreme 
points, stopping criteria for optimality at each iteration, 
the final interval of uncertainty in the line search and an 
upper limit on the maximum number of line searches con- 
ducted. Furthermore, the weights of the extreme points and 
the objective function evaluation are subject to roundoff 
errors. For a sensitive objective function, special care 
is necessary to insure convergence. Fortunately, we will 


confirm that the barrier function decomposition is a robust 
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procedure that does not require an optimal solution to the 
master problem at each iteration. Good results are obtained 
over a relative broad range of parameters. 

The comparison of the barrier algorithm RSD(B) with the 
penalty algorithm RSD(P) of Staniec [Ref. 4] reveals that 
they are very similar in their sequential structure. This 
similarity permits embedding both algorithms in the same 
computer program and creates the potential for devising a 
hybrid algorithm which takes advantage of each. The 
relationship between RSD(B) and RSD(P) is easily establish- 
ed. The dual estimates obtained from the penalty approach 


take, for some vector x and joint capacity constraint j, the 


form 

A -1,T 

E E (67) 
versus 

a,B = ní(o, - ax)Udji, (pb, -Ax)¿>0,j€J (68) 
for the barrier approach. However, the gradient of the 


penalty function at some vector x takes the same form as for 
the barrier function, namely 


VQ (h, xX) 


c + pia (69) 
versus 

A O 
Thus, besides the fact that the barrier function decomposi- 


tion requires an initial interior point, both methods use 


BER 


the same subproblems and master problem with different input 
for the dual estimates and different function evaluation 


routines in the line search. 
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IV. COMPUTATIONAL RESULTS 


In order to assess the capabilities of the algorithm 
RSD(B), we solve different versions of the test problem 
described in Chapter I. This problem suite has been exten- 
sively studied by Staniec [Ref. 4], using different al- 
gorithms. The optimal solution for four and ten commodity 
problems is available for comparisons, obtained by solving 
the problem (P) using the X-system [Ref. 32]. The four 
commodity problem (4H) has approximately 13,200 constraints, 
41,600 variables and optimal solution value 130,739,585. The 
ten commodity problem (10H) has about 33,000 constraints, 
104,000 variables and optimal solution value 169,532,339. 
For purposes of direct comparisons, the penalty algorithm 
RSD(P) of Staniec [Ref. 4] has been converted into our data 
structure and uses our computer program framework. Row Gr’) 
has been improved by taking advantage of the conjugate 
gradient modifications in the master problem. 

First, we will evaluate the performance of RSD(B) with 
different initial parameters in solving problem 4H. Subse- 
quently, the comparisons between RSD(P) and RSD(B) are 
presented for problems 4H and 10H. Possible modifications of 
algorithm RSD(B) are then discussed. Finally, we test with 
a 100 commodity problem having more than one million vari- 


ables and 300,000 constraints. 
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A. PERFORMANCE OF THE ALGORITHM RSD (B) 

The following results were obtained on an IBM 3033 AP 
under VM/CMS. The Central Processor Unit (CPU) utilization 
time is used as a performance measurement. An "e-optimality 
gap" is computed from the current upper and lower bounds as 
AAA cx or estimated as (V-V) /V 

Since the barrier function decomposition requires an 
interior starting point, we will first investigate the 
response of RSD(B) to different starting points. Let ca 


denote the maximum cost value in the network. While fixing 


n%=2*. c a = 0.5, and r = 7 in problem 4H, the choice 


max” 
of the parameter b establishing y" = b * y resulted in the 


optimality gaps listed in Table 1: 


Gap (5%) 


initial gap 66.76 51.93 39.42 2695 
iteration 7 WO), Zl 9787 SO 10.49 
iteration 16 4.60 3798 4.60 5.06 
iteration 24 Soo Ze Ss SUS IS 
iteration 32 2.03 1.42 1853 1586 
iteration 40 IS 0.99 eco 158 
iteration 48 0:78 0.68 0.89 0.89 
reduction b 0.80 0.85 90 0.95 


Table 1 : Response to Different Starting Points 


Problem 4H 
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If the parameter b is selected too low, the dual es- 
timates associated with the resulting interior starting 
point may not generate good, initial extreme points. If the 
starting value is chosen too high, the interior point gets 
too close to the boundary of the constraint region. The sub- 

sequent line search in the master problem is confined to a 


smaller search space of the constraint region resulting in 


slower convergence. Both values, b = 0.85 and b = 0.9, work 
well in the test problem. Further analysis will be based on 
b = 0.9. Good results with RSD were obtained for any value 


of r between 6 and 8. 

The lower bounds obtained with the algorithm RSD(B) are 
sensitive to the initial parameter hô., Recall that the dual 
estimates are computed as a- h(Ax-b,)~*. If the barrier 
parameter h? is relatively small, we approach the boundary 
very soon. As the slack on some jointly capacitated arcs 
almost vanishes, its reciprocal may take huge values result- 
ing in extreme dual estimates. This situation leads to large 
oscillations in subsequent solutions of the subproblem. The 
phenomenum is demonstrated in Figure 3 where hŷ = 0.5 * Se 
As h® is increased, this effect seems to disappear, as shown 
in Figures 4 through 6. It is interesting to observe that 
some good lower bounds are obtained under all conditions. 


The extreme result obtained for the small value of h? 


| 0, | | 
suggests that h” 2 Cmax 15 the better choice. 


41 


OPTIMUM 


0 100 250 300 


TIME (SECONDS) 


LOWER BOUND 
-5x107 0 5x10’ 1108 


—1x109 


Figure 3: V (fy) for h? = 0.5*Cmax, Problem 4H 
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Figure 4: V(p,) for h? = TC ae Problem 4H 
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Figure 5: V (Hy) for n? = 1.5*Cnm2x, Problem 4H 
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Figure 6: V (fy) for h? = 2*Cmax, Problem 4H 
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If we want to select an initial value h? that balances 
the linear part cx and the barrier term in the objective 
function, we have to take into account that |J,,| is not fixed 
but initially increasing. We do not have sufficient inform- 
ation about the size of the barrier term unless at least one 
pilot run has been conducted. Using the information that 
about 6% lof the arcs are finally containedAn secc O 


problem 4H, we would obtain a value of h < O II which 


max? 


is not the best choice, but may serve as a lower bound on hô. 
The differences obtained for the upper bounds are less 
significant. Figure 7 shows the sequence of values cx* 


obtained for the same choices of nO as before. 


1.8108 


1.8x 10% 


UPPER BOUNDS 


1.4108 


OPTIMUM 





1.2x10% 


0 100 220 300 
TIME (SECONDS) 


Figure 7: Response of CcxX to h, Problem 4H. 
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Any decrease in h creates a disturbance in the dual es- 
timates. Since the master problem is not solved optimally 
for each value of h, the algorithm may not always provide 
better lower bounds before the next update of h occurs. We 
are still able to get good lower bounds when a relatively 
moderate rate of decrease a = 0.5 is used. 

A solution trajectory for problem 4H is given in Figure 8 


using ho = 1.5 > c where a 0.6% optimal solution is 


max’ 


obtained within 200 seconds. Instead of plotting the 


eneu pper and lower bounds, the current values of V(t) 


k 


and cx” are shown. We observe the almost strictly decreasing 


ok 
sequence of cx 


although the master problem is not solved 
optimally. The improved upper bounds obtained by the 
Capacity allocation mechanism are denoted as cx’ and indi- 
cate that significant gains are obtained only at the initial 


stage. The differences between ox* 


and cx’ nearly vanish 
after about 100 seconds although slightly tighter upper 


bounds are produced throughout. 
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Figure 8: Trajectory of RSD(B), Problem 4H 


It seems to be typical for the barrier function that good 
upper bounds are obtained at an early stage, whereas the 
lower bounds trail behind. 

The trajectory of the objective value B(h,xX) together 
with its linear part cx“ is given in Figure 9. We find that 
B(h,x*) still approaches the optimum ex from below, the 
barrier term takes negative values over the whole range. The 
steps in its trajectory are due to the decreases of h by a 


factor of 075: 
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Figure 9: Values of the Objective Function, Problem 4H 


A direct comparison between the algorithms RSD(P) and 
RSD(B) shows that the barrier function decomposition is less 
effective in the solution of the smaller problem 4H but is 
competitive in the problem 10H. The solutions to problem 4H 
are compared in Figure 10 were ho = 0.001*c,7y for RSD(P) and 
for RSDB). 


Js 
h“ = E 
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Figure 10: Comparison of RSD(P) to RSD(B), Problem 4H 


If we use the same initial parameters for the solution of 
problem 10H, we obtain an interesting result. The trajec- 
tories are given in Figure 11. Obviously, the initial 
barrier and penalty parameters, respectively, are too low in 
both cases, yielding poor and oscillating lower bounds for 
RSD(B) versus bad upper bounds for RSD(P) due to insufficient 
penalty on the capacity violations. (This situation suggests 
investigation of a hybrid algorithm, incorporating both 


RSD(B) and RSD(P)). 
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Figure 11: RSD(P) Versus RSD(B) , Problem 10H, 


Same Initial Parameters 


The initial parameter hŷ has been increased until good 
results are obtained for both methods. The result is shown 
in Figure 12. We observe that the final values of v(pX) 
decay for both algorithms. We do not generate improving 
lower bounds. Apparently, the parameter h is updated to 
rapidly before the master problem generates good dual 
estimates. Staniec (Ref. 4] reported similiar poor lower 


bounding for problem 10H. 
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Figure 12: RSD(P) Versus RSD(B), Problem 10H, 


Best Parameters 


In summary, we find that RSD(B) provides good upper 
bounds at each iteration and does not depend on the capacity 
allocation routine as the penalty algorithm does. On: the 
other hand, the penalty algorithm provides better initial 
dual estimates, resulting in better lower bounds. Both 


algorithms converge to a good solution within a reasonable 


amount of time. 
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Bs POSSIBLE MODIFICATIONS OF THE ALGORITHM 

A first modification could be to postpone a further 
decrease of the barrier function parameter h until at least 
One better lower bound has been obtained compared to the 
bound generated for the previous value of h. Note, however, 
that the sequence of upper bounds is almost strictly decreas- 
ing from iteration to iteration ¿and controlled by the 
magnitude of h. A higher barrier parameter results in higher 
upper bounds. Therefore, this choice should depend on the 
problem at hand. If a fast approximation of the solution is 
desired, h should be decreased more rapidly. If a higher 
resolution is required and sufficient computer resources are 
available, a better solution of the master problem is 
necessary, yielding better intermediate lower bounds. 

Besides such a change of input parameters, a structural 
modification was also investigated. Since each capacity 
reallocation routine creates a subproblem solution that is 
feasible and at least as good as the current upper bound, 
this information can be passed to the master problem by using 
this solution as an additional "extreme" point. Experimenta- 
tion on the test problem showed however that no significant 
improvements were obtained. This may be because the dif- 
ference between the interior point solution of the master 


problem and the solution of the capacity reallocation 


gl 


procedure vanishes while the lower bounds are still trailing 
behind. Experimentation on other problems may yield dif- 
ferent results. 

The capacity allocation mechanism can be modified 
further. For RSD(B), capacity allocation for an interior 
point amounts to a redistribution of the available slack. As 
stated in equation (24), each commodity receives the same 


additional amount. A proportional allocation would be given 


Ko. + X03 (by - Ax) / (AX); 


ON Ypj — *p3 pj J J 


(71) 


The disadvantage of this procedure is that a commodity with 
zero flow on that arc gets zero capacity allocated. To 
overcome this drawback, a convex combination of both methods 
is possible: 

Zo3 (by 7 AX); / (AX) ; 
+ B, (by - AX) y E (72) 


Ypj T Epj + 8; 
where Bar Bo > 0 and Ba + Ba = 1. 

Experimentation with the test problem 4H showed improve- 
ments in the upper bounds in both cases, but since the lower 
bounds are still weak, the overall gain is not significant 
for this problem. 

RSD (B) would be superior to RSD(P) if we could establish 
better lower bounds. The dual estimates u are a function of 
the barrier parameter h and the slack on the joint capacity 


arcs. Different dual estimates are obtained if we use a 
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different value of h for each individual joint capacity arc. 
This approach leads to the weighted barrier function like the 
one proposed by Megiddo [Ref. 14] in general linear programm- 
OE 


B(h,x,w) = cx -~ h 2, AO A 


J yo 
where w is any real vector with positive components. Some 


(73) 


experiments were done with different weights. All arcs that 
are violated in the initial solution with completely relaxed 
joint capacity constraints are more likely to be tight in the 
opa mal solution: Therefore they are assigned a reduced 
weight to enable a faster approach to the boundary. Another 
weight factor used was proportional to the remaining slack on 
the corresponding arc. Unfortunately, neither attempt 
provided better solutions. 

The final modification is a hybrid barrier-penalty 
algorithm. Starting with the barrier function decomposition 
in problem 10H, we shift to the penalty algorithm after 16 
iterations using the extreme points generated so far. The 
main problem is the adjustment of the parameter h. Since we 
obtain good lower bounds for relatively small values of h, 
we reset h to the same initial value that we used in the 
independent approach. The results displayed in Figure 13 are 


unimpressive but further experimentation may improve results. 
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Figure 13 : Trajectory of the Hybrid Method 


Finally, we obtain a solution to the 100 commodity 
network flow problem 100H which is about 3.5% optimal within 
2650 seconds after 25 iterations and about 2.5% optimal 
within 3800 seconds after 38 iterations. The initial value 
hð has been increased to nO = 2ean Staniec [Ref. 4] 
reports a solution to this problem obtained by a penalty 
method that achieved 4 % optimality after 1000 seconds and 
finished with 1.5% optimality after 3000 seconds. Our method 


did not show the same performance. The solution trajectory 


as shown in Figure 14 seems to indicate that the upper bounds 


54 


are already very tight whereas the lower bounds are again 
poor. The objective function value (not shown) is still less 
than the lower bound, a further decrease of h is necesarry 


for convergence. 
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Figure 14 : Solution Trajectory, Problem 100H 


The distribution of the CPU time between the master 
problem and the subproblems changes as we increase the number 
of commodities. Although the solution of each pure network 
flow problem requires less than one second CPU time per 


commodity in the average for all test problems, we find that 


the subproblems are expensive to solve and consume the 
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largest portion of the total CPU time TAC eo neers 
shown in Table 2. The amount of time not used by the master 
problem or the subproblem is mostly consumed in the 


input/output operations. 


Problem 4H Problem 10H Problem 100H 


Master Problem 25.5% 11.2% 7.1% 
Subproblem 49.9% 66.9% 86.3% 
Other 24.6% 21.9% 6.6% 


Table 2: Distribution of CPU Time for Test 


Problems of Increasing Size 


However, the subproblem solves for each commodity indepen- 
dently and more than one pure network flow problem could be 
solved at a time. This feature makes our method highly 
parallel with a potential reduction of nearly 1/|P| in 


elapsed computation time for the subproblem. 


56 


V. CONCLUSIONS 


The method of logarithmic barrier function decomposition 
is a useful tool for the solution of large-scale multicom- 
modity flow problems. The algorithm RSD(B) is a variant of 
the price-directive decomposition method. It generates a 
sequence of interior points which provide intermediate 
feasible solutions while converging towards the optimum. The 
technique of restricted simplicial decomposition (RSD) proves 
to be useful in the solution of the nonlinear master problem. 
However, RSD(B) seems to be robust and does not require an 
optimal solution to the master problem. It is competitive 
with penalty decomposition methods and relatively easy to 
implement. Since it achieves good feasible solutions at an 
early stage, it may be used as a starting technique in other 
algorithms like the hybrid barrier-penalty technique. The 
use of RSD(B) in other large-scale multicommodity flow 


problems is recommended. 
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